function [theta0,theta0TransReduced,setup]   = loadThetaStart(fileStart,setup)

if exist([fileStart,'.txt'],'file')
    disp(['Loading solution:  Approximation order = ',num2str(setup.appOrder)])
    % Dimension of theta: ny * numDim2
    theta0 = load([fileStart,'.txt']);
else    
    % Starting values from the perturbation approximation
    modelPer = setup.modelPer;
    nx       = modelPer.nx;
    if setup.appOrder == 1
        theta0 = [modelPer.g0+0.5*modelPer.gss   modelPer.gx.*repmat(setup.scaleX_in_g',size(modelPer.gx,1),1)];
    elseif setup.appOrder == 2
        theta0 = zeros(modelPer.ny,1+nx+nx*(nx+1)/2);
        theta0(:,1:1+nx) = [modelPer.g0+0.5*modelPer.gss  modelPer.gx.*repmat(setup.scaleX_in_g',size(modelPer.gx,1),1)];
        idx = 1+nx;
        for i=1:nx
            for j=i:nx
                idx = idx + 1;
                theta0(:,idx) = 0.5*modelPer.gxx(:,i,j)*(setup.scaleX_in_g(i,1)*setup.scaleX_in_g(j,1));
            end
        end
    elseif setup.appOrder == 3
        theta0 = zeros(modelPer.ny,1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6);
        theta0(:,1:1+nx) = [modelPer.g0+0.5*modelPer.gss  (modelPer.gx+3/6*modelPer.gssx).*repmat(setup.scaleX_in_g',size(modelPer.gx,1),1)];
        idx = 1+nx;
        for i=1:nx
            for j=i:nx
                idx = idx + 1;
                theta0(:,idx) = 0.5*modelPer.gxx(:,i,j)*(setup.scaleX_in_g(i,1)*setup.scaleX_in_g(j,1));
            end
        end
        for alfa1=1:nx
            for alfa2=alfa1:nx
                for alfa3=alfa2:nx
                    idx = idx + 1;
                    theta0(:,idx) = 1/6*modelPer.gxxx(:,alfa1,alfa2,alfa3)*...
                        (setup.scaleX_in_g(alfa1,1)*setup.scaleX_in_g(alfa2,1)*setup.scaleX_in_g(alfa3,1));
                end
            end
        end
    end
end

% Setting the approximation order per variable
modelPer    = setup.modelPer;
nx          = modelPer.nx;
appOrderVar = setup.appOrderVar;
select      = ones(size(theta0));
for i=1:size(select,1)
    if appOrderVar(i) == 1
        select(i,1+nx+1:end) = 0;
    elseif appOrderVar(i) == 2
        select(i,1+nx+nx*(nx+1)/2+1:end) = 0;
    end
end
theta0Trans        = theta0';
theta0TransReduced = theta0Trans(select'==1);
setup.select = select;
end


